# ============================================================================
# Integrated Soil Seismic Response Analysis - Simplified Version
# OpenSees Tcl Script for 2D Soil Model Dynamic Analysis
# ============================================================================

# Clear existing model
wipe

puts "Integrated Soil Seismic Analysis Started"

# ============================================================================
# 1. Configuration Parameters
# ============================================================================

# Earthquake file settings
set earthquakeFile "GroundMotions/earthquake.txt"
set timeStep 0.001             ;# Analysis time step (sec) - must match earthquake data
set analysisDuration 15.0      ;# Total analysis duration (sec) - test shorter duration

# Results output directory
set timestamp [clock format [clock seconds] -format "%Y%m%d_%H%M%S"]
set outputDir "Results/$timestamp"
file mkdir $outputDir
file mkdir "$outputDir/Acceleration"
file mkdir "$outputDir/Displacement"

# ============================================================================
# 2. Finite Element Model Creation
# ============================================================================

model basic -ndm 2 -ndf 2

# Geometry parameters
set soilWidth    20.0    ;# Soil width (m)
set soilHeight   10.0    ;# Soil height (m)
set nEleX        20      ;# Number of elements in X direction
set nEleY        10      ;# Number of elements in Y direction

# Calculate mesh spacing
set dx [expr $soilWidth/$nEleX]
set dy [expr $soilHeight/$nEleY]

# Material parameters
set rho        1800.0   ;# Soil density (kg/m3)
set E          50.0e6   ;# Elastic modulus (Pa)
set nu         0.3      ;# Poisson ratio

# Create nodes
set nodeTag 1
for {set j 0} {$j <= $nEleY} {incr j} {
    for {set i 0} {$i <= $nEleX} {incr i} {
        set x [expr $i * $dx]
        set y [expr $j * $dy]
        node $nodeTag $x $y
        incr nodeTag
    }
}

# Define material - Linear Elastic
# For plane strain condition: ElasticIsotropic material
nDMaterial ElasticIsotropic 1 $E $nu $rho

# Create elements
set eleTag 1
for {set j 0} {$j < $nEleY} {incr j} {
    for {set i 0} {$i < $nEleX} {incr i} {
        set n1 [expr $j*($nEleX+1) + $i + 1]
        set n2 [expr $j*($nEleX+1) + $i + 2]  
        set n3 [expr ($j+1)*($nEleX+1) + $i + 2]
        set n4 [expr ($j+1)*($nEleX+1) + $i + 1]
        
        element quad $eleTag $n1 $n2 $n3 $n4 1.0 "PlaneStrain" 1
        incr eleTag
    }
}

# Boundary conditions - Based on modelEN.tcl reference
# Bottom: fully fixed (bedrock) - complete fixity like modelEN nodes 1 and 2
for {set i 1} {$i <= [expr $nEleX+1]} {incr i} {
    fix $i 1 1  ;# X and Y translation fixed (bedrock condition)
}

# Left and Right sides: Fixed in horizontal direction to prevent sliding
# This follows the modelEN approach of constraining base nodes
for {set j 1} {$j <= $nEleY} {incr j} {
    # Left side nodes
    set leftNode [expr $j*($nEleX+1) + 1]
    fix $leftNode 1 0  ;# X-direction fixed, Y free
    
    # Right side nodes  
    set rightNode [expr $j*($nEleX+1) + ($nEleX+1)]
    fix $rightNode 1 0  ;# X-direction fixed, Y free
}

# ============================================================================
# 3. Gravity Equilibrium Analysis
# ============================================================================

puts "Starting gravity equilibrium analysis..."

# Gravity load pattern
pattern Plain 1 Linear {
    for {set j 0} {$j <= $nEleY} {incr j} {
        for {set i 0} {$i <= $nEleX} {incr i} {
            set nodeNum [expr $j*($nEleX+1) + $i + 1]
            
            if {$j == 0} continue
            
            # Calculate tributary area
            set tributaryArea 0.0
            if {($i == 0 || $i == $nEleX) && ($j == $nEleY)} {
                set tributaryArea [expr $dx * $dy * 0.25]
            } elseif {($i == 0 || $i == $nEleX) && ($j > 0 && $j < $nEleY)} {
                set tributaryArea [expr $dx * $dy * 0.5]
            } elseif {($i > 0 && $i < $nEleX) && ($j == $nEleY)} {
                set tributaryArea [expr $dx * $dy * 0.5]
            } else {
                set tributaryArea [expr $dx * $dy]
            }
            
            set bodyForce [expr -$rho * 9.81 * $tributaryArea]
            if {$tributaryArea > 0.0} {
                load $nodeNum 0.0 $bodyForce
            }
        }
    }
}

# Gravity analysis settings
constraints Plain
numberer Plain
system BandGeneral
test NormDispIncr 1.0e-6 25
algorithm Newton
integrator LoadControl 0.05
analysis Static

# Perform gravity analysis
set ok [analyze 20]
if {$ok != 0} {
    integrator LoadControl 0.02
    set ok [analyze 50]
}

if {$ok == 0} {
    puts "Gravity equilibrium analysis completed"
} else {
    puts "Warning: Gravity analysis not fully converged"
}

loadConst -time 0.0

# Clear analysis object before setting up dynamic analysis
wipeAnalysis

# ============================================================================
# 4. Setup Recorders - X direction only, no time column
# ============================================================================

set topCenterNode [expr $nEleY*($nEleX+1) + ($nEleX+1)/2 + 1]

# Main recorders - X direction only, no time column
recorder Node -file "$outputDir/Displacement/TopDisplacement.out" -node $topCenterNode -dof 1 disp
recorder Node -file "$outputDir/Acceleration/TopAcceleration.out" -node $topCenterNode -dof 1 accel

# Top surface displacement - X direction only
set topNodes {}
for {set i 0} {$i <= $nEleX} {incr i} {
    set nodeNum [expr $nEleY*($nEleX+1) + $i + 1]
    lappend topNodes $nodeNum
}
recorder Node -file "$outputDir/Displacement/TopSurfaceDisplacement.out" -node {*}$topNodes -dof 1 disp
recorder Node -file "$outputDir/Acceleration/TopSurfaceAcceleration.out" -node {*}$topNodes -dof 1 accel

# ============================================================================
# 5. Seismic Analysis
# ============================================================================

puts "Starting seismic analysis..."

# Check earthquake file
if {![file exists $earthquakeFile]} {
    puts "Error: Earthquake file not found: $earthquakeFile"
    exit
}

# Read earthquake data (keep original values in g)
set inFile [open $earthquakeFile r]
set earthquakeData {}
set timePoints {}
set currentTime 0.0
set dataTimeStep $timeStep

while {[gets $inFile line] >= 0} {
    set line [string trim $line]
    # Skip empty lines and comments
    if {$line == "" || [string index $line 0] == "#"} {
        continue
    }
    
    # Clean up the line - remove extra spaces and tabs
    # regsub -all {\s+} $line " " line
    set data [split $line]
    
    if {[llength $data] >= 2} {
        # Format: time acceleration
        set timeStr [lindex $data 0]
        set accelStr [lindex $data 1]
        
        # Check if values are numeric - keep acceleration in g units
        if {[string is double -strict $timeStr] && [string is double -strict $accelStr]} {
            set time $timeStr
            set accel $accelStr  ;# Keep in g units
            
            lappend timePoints $time
            lappend earthquakeData $accel
        } else {
            puts "Warning: Skipping non-numeric data: $line"
        }
    } elseif {[llength $data] == 1} {
        # Format: acceleration only (assume uniform time step)
        set accelStr [lindex $data 0]
        
        # Check if value is numeric - keep acceleration in g units
        if {[string is double -strict $accelStr]} {
            set accel $accelStr  ;# Keep in g units
            
            lappend timePoints $currentTime
            lappend earthquakeData $accel
            set currentTime [expr $currentTime + $dataTimeStep]
        } else {
            puts "Warning: Skipping non-numeric data: $line"
        }
    }
}
close $inFile

# Check if we have valid data
if {[llength $earthquakeData] == 0} {
    puts "Error: No valid earthquake data found in file: $earthquakeFile"
    puts "Please check the file format. Expected formats:"
    puts "  1) Two columns: time acceleration"  
    puts "  2) Single column: acceleration (uniform time step assumed)"
    exit
}

puts "Loaded [llength $earthquakeData] earthquake data points"
puts "Time range: [lindex $timePoints 0] to [lindex $timePoints end] seconds"
puts "Max acceleration: [lindex [lsort -real $earthquakeData] end] g (will be multiplied by 9.8)"

# Calculate data time step
if {[llength $timePoints] >= 2} {
    set dt [expr [lindex $timePoints 1] - [lindex $timePoints 0]]
} else {
    set dt $timeStep
}

puts "Analysis time step: $timeStep sec, Data time step: $dt sec"

# Create earthquake load using the correct method (like the reference)
# Write acceleration data in g units to temporary file
set tempFile "$outputDir/temp_earthquake.dat"
set outFile [open $tempFile w]
foreach accel $earthquakeData {
    puts $outFile $accel
}
close $outFile

# Create time series with factor to convert g to m/s²
timeSeries Path 1 -dt $dt -filePath $tempFile -factor 9.8
pattern UniformExcitation 2 1 -accel 1

# Setup damping
set dampRatio 0.05
set omega1 [expr 2*3.14159*1.0]   ;# 1Hz
set omega2 [expr 2*3.14159*10.0]  ;# 10Hz
set a0 [expr $dampRatio*2*$omega1*$omega2/($omega1+$omega2)]
set a1 [expr $dampRatio*2/($omega1+$omega2)]
rayleigh $a0 0.0 0.0 $a1

# Dynamic analysis settings
constraints Plain
numberer Plain

# Test both solvers for performance comparison
puts "=== Testing BandGeneral solver ==="
set startTime [clock clicks -milliseconds]
system BandGeneral
test NormDispIncr 1.0e-6 25
algorithm Newton
integrator Newmark 0.5 0.25
analysis Transient

# Test run with BandGeneral (10 steps only for comparison)
set bandTime 0
for {set i 1} {$i <= 10} {incr i} {
    set testStart [clock clicks -milliseconds]
    set ok [analyze 1 $timeStep]
    set testEnd [clock clicks -milliseconds]
    set bandTime [expr $bandTime + ($testEnd - $testStart)]
    if {$ok != 0} break
}
puts "BandGeneral: 10 steps took ${bandTime} milliseconds"

# Clear and setup GPU solver
wipeAnalysis
puts "=== Testing GPUSolver ==="
system GPUSolver

test NormDispIncr 1.0e-6 25
algorithm Newton
integrator Newmark 0.5 0.25
analysis Transient

# Test run with GPUSolver (10 steps only for comparison)
set gpuTime 0
for {set i 1} {$i <= 10} {incr i} {
    set testStart [clock clicks -milliseconds]
    set ok [analyze 1 $timeStep]
    set testEnd [clock clicks -milliseconds]
    set gpuTime [expr $gpuTime + ($testEnd - $testStart)]
    if {$ok != 0} break
}
puts "GPUSolver: 10 steps took ${gpuTime} milliseconds"

# Performance comparison
if {$bandTime > 0 && $gpuTime > 0} {
    set speedup [expr double($bandTime) / $gpuTime]
    puts "Performance comparison: GPU/CPU speedup = [format %.2f $speedup]x"
    if {$speedup < 1.0} {
        puts "WARNING: GPU is slower than CPU for this problem size"
        puts "GPU overhead may dominate for small problems"
    }
}

# Continue with full analysis using the faster solver
wipeAnalysis
if {$bandTime < $gpuTime} {
    puts "Continuing full analysis with BandGeneral (faster)"
    system BandGeneral
} else {
    puts "Continuing full analysis with GPUSolver"
    system GPUSolver
}

test NormDispIncr 1.0e-6 25
algorithm Newton
integrator Newmark 0.5 0.25
analysis Transient

# Perform dynamic analysis
set numSteps [expr int($analysisDuration / $timeStep)]
set progressInterval [expr $numSteps / 20]  ;# Show progress every 5%

puts "Starting dynamic analysis with $numSteps steps..."

set ok 0
for {set step 1} {$step <= $numSteps} {incr step} {
    set ok [analyze 1 $timeStep]
    
    # Show progress every 5%
    if {$step % $progressInterval == 0} {
        set progress [expr int(100.0 * $step / $numSteps)]
        set currentTime [expr $step * $timeStep]
        puts "Progress: ${progress}% (${currentTime}s / ${analysisDuration}s)"
    }
    
    if {$ok != 0} {
        algorithm ModifiedNewton
        set ok [analyze 1 $timeStep]
        if {$ok != 0} {
            puts "Analysis failed at step $step (time: [expr $step * $timeStep]s)"
            break
        }
        algorithm Newton
    }
}

if {$ok == 0} {
    puts "Dynamic analysis completed"
} else {
    puts "Dynamic analysis partially completed"
}

# ============================================================================
# 6. Simple Results Analysis
# ============================================================================

# Read maximum displacement (X direction only, no time column)
if {[file exists "$outputDir/Displacement/TopDisplacement.out"]} {
    set dispFile [open "$outputDir/Displacement/TopDisplacement.out" r]
    set maxDispX 0.0
    
    while {[gets $dispFile line] >= 0} {
        set dispX [string trim $line]
        
        if {[string is double -strict $dispX]} {
            if {abs($dispX) > abs($maxDispX)} {set maxDispX $dispX}
        }
    }
    close $dispFile
    
    puts "Max X displacement: [format %.4f $maxDispX]m"
}

# Read maximum acceleration (X direction only, no time column)
if {[file exists "$outputDir/Acceleration/TopAcceleration.out"]} {
    set accelFile [open "$outputDir/Acceleration/TopAcceleration.out" r]
    set maxAccelX 0.0
    
    while {[gets $accelFile line] >= 0} {
        set accelX [string trim $line]
        
        if {[string is double -strict $accelX]} {
            if {abs($accelX) > abs($maxAccelX)} {set maxAccelX $accelX}
        }
    }
    close $accelFile
    
    puts "Max X acceleration: [format %.4f $maxAccelX]m/s²"
}

# Generate summary report
set summaryFile [open "$outputDir/Summary.txt" w]
puts $summaryFile "Soil Seismic Response Analysis Results"
puts $summaryFile "Model: ${nEleX}x${nEleY} elements, ${soilWidth}mx${soilHeight}m"
puts $summaryFile "Analysis steps: $step/$numSteps"
if {[info exists maxDispX]} {
    puts $summaryFile "Max X displacement: [format %.4f $maxDispX]m"
}
if {[info exists maxAccelX]} {
    puts $summaryFile "Max X acceleration: [format %.4f $maxAccelX]m/s²"
}
puts $summaryFile ""
puts $summaryFile "Output files:"
puts $summaryFile "- X displacement data: Displacement/ folder"
puts $summaryFile "- X acceleration data: Acceleration/ folder"
puts $summaryFile "- Data format: Single column (no time column)"
close $summaryFile

puts "Analysis completed, results saved in: $outputDir"